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Abstract 



We consider a standard splitting algorithm for the rare-event simulation of overflow 
probabilities in any subset of stations in a Jackson network at level n, starting at a fixed 
initial position. It was shown in [8j that a subsolution to the Isaacs equation guarantees 
that a subexponential number of function evaluations (in n) suffice to estimate such 
overflow probabilities within a given relative accuracy. Our analysis here shows that in 
fact O (n 2 ' 3+1 ) function evaluations suffice to achieve a given relative precision, where 
j3 is the number of bottleneck stations in the network. This is the first rigorous analysis 
that allows to favorably compare splitting against directly computing the overflow 
probability of interest, which can be evaluated by solving a linear system of equations 
£T| ■ with 0(n d ) variables. 

o 

q : 1 Introduction 

o' 

The development of rare-event simulation algorithms for overflow probabilities in stable open 
Jackson networks has been the subject of a substantial amount of papers in the literature 
during the last decades (see Section 2 for the specification of an open Jackson network). 
A couple of early references on the subject are [20] and [lj. Subsequent work which has 
also been very influential in the development of efficient algorithms for overflows of Jackson 
networks include [22], [32], [13], [33], [33], [9], [19], [TT] and [8]. The survey papers of [33] 
and [7] provide additional references on this topic. 

The two most popular approaches that are applied to the construction of efficient rare- 
event simulation algorithms are importance sampling and splitting (see [3]). Importance 
sampling involves simulating the system under consideration (in our case the Jackson net- 
work) according to a different set of probabilities in order to induce the occurrence of the 
rare event. Then, one attaches a weight to each simulation corresponding to the likelihood 
ratio of the observed outcome relative to the nominal / original distribution. In splitting, on 
the other hand, there is no attempt to bias the behavior of the system. Instead, the idea is 
to decompose the occurrence of the rare event of interest (in our case overflow in a Jackson 
network) into a sequence of nested "milestone" events whose subsequent occurrence is not 
rare. The rare event occurs when the last of the milestone events occurs. The idea is to 
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keep splitting the particles as they reach subsequent milestones. Of course, each particle is 
attached a weight corresponding to the total number of times it has split so that the overall 
estimation (which is the sum of the weights corresponding to the particles that make it to 
the last milestone) provides an unbiased estimator of the probability of interest. 

The most popular performance measure for efficiency analysis of rare-event simulation 
algorithms for Jackson networks corresponds to that of "asymptotic optimality" or "weak 
efficiency" . In order to both explain the computational complexity implied by this notion and 
to put in perspective our contributions let us discuss the class of problems we are interested 
in. Starting from any fixed state, we consider the problem of computing the probability that 
the total number of customers in any fixed set of stations in the network reaches level n prior 
to reaching the origin. That is, the probability that the sum of the queue lengths in any 
given set of stations reaches level n within a busy period. The number of stations in the 
whole network is assumed to be d and the number of bottleneck stations (i.e. stations with 
the maximum traffic intensity in equilibrium) is (3. 

Weak efficiency guarantees that a subexponential number of replications (as a function of 
the overflow level, say n) suffices to compute the underlying overflow probability of interest 
within a given relative accuracy. In contrast, as we shall explain in Section 2, overflow 
probabilities in the setting of Jackson networks can be computed by solving a linear system 
of equations with 0(n d jj unknowns. It is well known that Gaussian elimination then takes 
O (n M ) operations to find an exact solution to such linear system. Moreover, since in our 
case the associated linear system has some sparsity properties the linear equations can be 
solved in as many as O (n 3d ~ 2 ) operations (see the discussion in Section 2). Our analysis for 
the solution of the associated linear system of equations is not intended to be exhaustive. 
Our objective is simply to make the point that naive Monte Carlo (which indeed takes an 
exponential number of replication in n to achieve a given relative accuracy) is not the natural 
benchmark that one should be using in order to test the performance of an efficient simulation 
estimator for overflows in Jackson networks. Rather, a more natural benchmark is the 
application of a straightforward method for solving the associated system of linear equations. 
It would be interesting to provide a detailed study of various methods for solving linear system 
of equations (such as multigrid procedures) that are suitable for our environment and even 
combine them with the ideas behind efficient simulation procedures. This, however, would 
be the subject of an entire paper and therefore is left as a topic for future research. 

Our goal here is to analyze a class of splitting algorithms similar to those introduced in 
[22] for the evaluation of overflow probabilities at level n. Further analysis was given in [8], 
where the authors provide necessary and sufficient conditions for the design of the "milestone 
events" in order to achieve subexponential complexity in n. Our contribution is to show 
that if the milestone events are properly placed as suggested by pj, the splitting algorithm 
requires O (n 2 ^ +1 ) function evaluations to achieve a fixed relative error. Since clearly the 
number of bottleneck stations (3 is at most d, the complexity of splitting is O (n 2d+1 ), which 
is substantially smaller than that of the direct solution of the associated linear system. Our 
analysis therefore provides theoretical justification for the superior performance observed 



^^Given two non-negative functions /(■) and <?(■), we say / (n) = 0{g{n)) if there exists c, no 6 (0,oo) 
such that / (n) < eg (n) all n > uq. 



when applying splitting algorithms compared to directly solving the associated linear system. 

We believe that our results shed light into the type of performance that can be expected 
when applying particle algorithms beyond the setting of Jackson networks. This feature 
should be emphasized, specially given the fact that a linear time algorithm for computing 
overflows in Jackson networks has been developed very recently (see [1]). Contrary to particle 
methods, which are versatile and that can in principle be applied in great generality, the 
algorithm in [3] takes advantage of certain properties of Jackson networks which are not 
shared by all classes of systems. 

In addition, our results also provide interesting connections to recent performance analysis 
studied in the context of state- dependent importance sampling algorithms for a class of 
Jackson networks. These connections might eventually help guide the users of rare event 
simulation algorithms decide when to apply importance sampling or splitting. For instance, 
consider the overflow at level n of the total population of a tandem network with d stations. 
The work of [S] proposes an importance sampling estimator based on the subsolution of an 
associated Isaacs equation. In particular, [H] shows that if exponential tiltings are applied 
using the gradient of the associated subsolution as the tilting parameter (depending on the 
current state), the corresponding algorithm is weakly efficient. Turns out that there are many 
subsolutions that can be constructed varying certain so-called "mollification parameters" . A 
recent analysis based on Lyapunov inequalities given in [6] shows that a natural selection 
of mollification parameters guarantees O (n 2< - d_/3 * )+1 ) function evaluations to achieve a given 
relative error. Our analysis here therefore guarantees that one can achieve a running time of 
order O (n d+1 ) if one chooses importance sampling when there are more than d/2 bottleneck 
stations in the network and splitting if there are less than d/2 bottleneck stations. Although 
our analysis is still not sharp we believe that our results provide a significant step in order 
to understand the connections between splitting and importance sampling. 

The rest of the paper is organized as follows. A brief discussion on complexity and 
efficiency considerations is given in Section [2j Then we discuss the necessary large deviations 
asymptotics for Jackson networks required for our analysis in Section [3j The introduction 
of the splitting algorithm as well as connections to the theory developed in [8] is given in 
Section |U Our complexity analysis is finally given in Section 

2 Complexity and Efficiency 

We shall review concepts of efficiency and complexity in rare event simulation. We start 
our discussion in the context of a generic class of rare event simulation problems. Consider 
a sequence of events {E n ,n = 1,2...} with p n = P (E n ) — > as n — )■ oo (Without loss 
of generality, we might assume that p n — > exponentially fast as n /* 0.) The design of 
an efficient rare-event simulation algorithm is typically associated with the construction of 
an unbiased estimator, say p n , such that p n = E[p n ]. A number of m i.i.d. replications 
{pn i ■■■ipn } is produced, the average of which forms an estimate of p n , namely 

in 



p n {m) 

m 

3=1 



—J2p 

rn < • 



By virtue of Chebyshev's inequality we obtain the following property for the relative error, 
\p n (m) — p n \/Pn, of the estimate 

Var(j3 n ) 



P(|p„(m) -Pn\/Pn > e) < 



mp^e 2 



Hence, for a pre-determined upper bound e of relative error, if we choose the number of 
replications m such that 

m > e- 2 5- 1 (cv n ) 2 , (2) 

where cv 2 = Var (p n ) jp\ is the squared coefficient of variation of p n , we can guarantee that 
the relative error is no larger than e with probability at least 1 — 5. 

Equation ([2]) stipulates that m needs to grow at least at the same rate as cv 2 does in 
order to keep the relative error within a desirable threshold. If cv 2 grows at a subexponen- 
tial rate (i.e. if log (cv n ) = o(\ogp n ) , asn / oo) the estimator is said to be asymptotic 
optimality, logarithmical efficient or weakly efficient. In this case, the number of replication 
needs to increase subexponentially in n to achieve a prescribed level of relative accuracy. 
The name "asymptotic optimality" is derived from the fact that weak efficiency implies that 
the exponential rate of decay to zero of the Ep^ coincides with that of p 2 n and therefore is 
maximal (by virtue of Jensen's inequality). 

Obviously, one has to keep in mind that weak efficiency measures the optimality of the 
estimator for a given level of computational budget. Coming to the splitting algorithm, it 
is apparent that the computational effort varies drastically with the degree of splitting per- 
formed; one must therefore take into account the cost involved in generating each replication 
of p n . We measure such cost in terms of the number of elementary function evaluation which 
we will take to be simple addition, multiplication, comparison, and the generation of a single 
uniform random variable. When we incorporate the computational cost per replication of 
the estimator, <^j says that the total number of function evaluations needed has to keep pace 
with the work-normalized squared coefficient of variation, i.e., cv 2 , ■ N n , where N n is the cost 
per replication of p n . We will show in section [5] that N n is closely related to the expected 
total number of the survival particles in a single run of the Splitting algorithm. 

Coming back to the setting of Jackson networks. It is important to recognize that over- 
flow probabilities in such a setting can be obtained by solving a system of linear equations. 
So, a reasonable benchmark for any simulation based algorithm to be regarded "efficient" is 
indeed how fast one can solve such a system of linear equation by a direct procedure. Jackson 
networks are basically multidimensional simple random walks with constrained behavior on 
the boundaries. In particular, they are Markov chains living on a countable state-space. The 
overflow probabilities can be conveniently expressed as first passage time probabilities, which 
in turn can be characterized as the solution to certain linear system of equations thanks to its 
countable state-space Markov chain structure. We shall quickly review how to obtain such 
linear system for a generic Markov chain Q = {Qk : k > 0} living on a countable state-space 
S with transition matrix {K (x, y) : x,y E S}. Let A, B be two disjoint subsets of S, define 
T A = mi{k > : X E A}, T B = ml{k > : X E B} and put p(x) = F x (T A < T B ). A 
simple conditioning argument on the first transition leads to 

p( x ) = ^2 K ( x ,y)p(y) ( 3 ) 

yes 



subject to the boundary conditions 

p (x) = 1 for x G A, p(x) = for x G B. 

In fact, p(-) is the minimum non-negative solution to the above system (see [5]). 

Now, if Q describes the state of the embedded discrete time Markov chain corresponding 

to a Jackson network with d stations then S = Z+. The transition dynamics of a Jackson 

network are specified as follows (see [21] p. 92). Inter-arrival times and the service times 

are all independent and exponentially distributed random variables. The arrival rates are 

given by the vector A = (Ai,...,Ad) and service rates are given by fi = (/^i, ...,fid) ■ (By 

convention all of the vectors in this paper are taken to be column vectors and T denotes 

transposition.) A job that leaves station i joins station j with probability Pij and it leaves 

the system with probability 

d 

Pifi = 1 - / „ Pj,y 

The matrix P = {Pi.j : 1 < i,j ' < d} is called the routing matrix. We shall consider open 
Jackson networks, which satisfy the following conditions: 

i) Vi, either A» > or Xj 1 Pj 1 j 2 ...Pj k i > for some ji, ...,jk- 
ii) Vi, either P i0 > or P i j 1 Pj 1 j 2 ...Pj k0 > for some ji, ...,jk- 
iii) The network is stable (i.e. a stationary distribution exists). 

These conditions simply require that each station will receive jobs either directly from the 
outside or routed from other stations, and each job will leave the system eventually. Our main 
interest lies in the evaluation of p n (x) assuming that B = {0} and A n = {y : v T y = n} where 
v is a binary vector which encodes a particular subset of the network (i.e., the z-th position of 
the vector v is 1 if station i falls in the subset of interest, and otherwise). We shall denote 
by V (x) = x T v the mapping recording the total population in the stations corresponding to 
the vector v. The case in which v — 1 = (1, 1, ..., 1) corresponds to the total population of 
the system. So, p n (x), or more precisely p\ (x), corresponds to the overflow probability in 
the subset encoded by v within a busy period and starting from x and that. In this setting, 
it follows (as we shall review in the next section) that p\ (x) — > exponentially fast in n as 
n /*• oo and the system of equations ([3D has O (n d ) unknowns. Gaussian elimination requires 
O (n 3d ) function evaluations to find the solution of such system. But since each state of the 
Markov chain in this case has possible interactions with only a small fraction of the entire 
state-space, it is therefore possible to permutate the states (say in lexicographic order) so 
that the system is banded (i.e. the associated matrix is sparse in the sense that its non-zero 
entries fall to a diagonal band.) One can show that the bandwidth is O (n d_1 ) , and therefore 

solving such a banded linear system requires O (n d ■ (n d_1 ) J = O (n M ~ 2 ) operations (see, 
e.g., 0). 



Estimators that possess weak efficiency (in a work-normalized sense) are guaranteed to 
run at subexponential complexity. When comparing to the above polynomial algorithms 
of solving systems of linear equations, the efficiency analysis of such estimators appears 
to be insufficient. We will show in later analysis that the multilevel Splitting algorithm 
suggested by Dean and Dupuis jHj, applied to estimate the overflow probabilities in Jackson 
networks, requires fewer function evaluations than directly solving the associated system of 
linear equations. 

3 Jackson Networks: Notation and Properties 

As we mentioned in the previous section, a Jackson network is encoded by two vectors of 
arrival and service rates, A = (Ai, ■■■,Xd) and p = (p\, ...,pd) , together with a routing 
matrix P = {Py : i <i,j < d}. Without loss of generality, we assume that ^2 i=1 (Aj + Pi) = 
1. The network is assumed to be open and stable so conditions i), ii), and iii) described in 
the previous section are in place. 

Given the stability assumption, the system of equations given by 

d 

<f>i = K + ^2<f>jPji, Vz = l,2,...,d (4) 

i=i 

admits a unique solution = X T (I — P) (see [3]). The traffic intensity at station i in the 
system in equilibrium is given by pi which is defined by 

Hi Hi 

and satisfies pi G (0,1) for all i = 1,2, ...,d. Define p* = maxi<i<dPi and let (3 be the 
cardinality of the set {i : pi = p*}. 

We shall study the queueing network by means of the embedded discrete time Markov 
chain Q = {Q(k) : k > 0}, where Q(k) = (Qi(k), . . . ,Qd(k)). For each k, Qi(k) represents 
the number of customers in station % immediately after the kth transition epoch of the system. 
As mentioned before, the process Q lives in the space S = Zi, 

Let V (y) = y T v be the total population in the stations corresponding to the binary vector 
v. We are interested in the overflow probability in any given subset of the Jackson network. 
More precisely, we wish to estimate 

p^ = P { total population in stations encoded by v reaches 
n before returning to 0, starting from 0}. 

In turn, p^ can be expressed in terms of the following stopping times, 

T {0} ^ inf{£; > 1 : Q (k) = 0}, 
T^M{k> 1: V(Q(k)) >n}. 



Indeed, if we use the notation P x (-) = P(-|Q(0) = x) then we can rewrite p\ as 

pl=F (TV<T {0} ). (6) 

Similarly, 

pl(x)=F x (TV<T {0} ). (7) 

The asymptotic analysis of p^ (x) can be studied by means of large deviations theory. 
We shall indicate how this theory can be applied to specify an efficient splitting algorithm 
in the next section. In the mean time, let us provide a representation for the dynamics of 
the queue length process that will be convenient in order to motivate the elements of the 
efficient splitting algorithm that we shall analyze. 

As we mentioned earlier, Jackson networks are basically constrained random walks. The 
constraints arise because the number of customers in each station must be non-negative. 
Thinking about Jackson networks as constrained random walks facilitates the introduction 
and motivation of the necessary large deviations elements behind the description of the 
splitting algorithm. In order to specify the dynamics of the embedded discrete time Markov 
chain in terms of a random walk type representation we need to introduce notation which 
will be useful to specify the transitions at the boundaries induced by the non-negativity 
constraints. 

The state-space Z+ can be partitioned in 2 d different regions which are indexed by all 
the subsets E C {1, . . . , d}. The region encoded by a given subset E is defined as 

d E = {z e Z d + : Zi = 0, i e E, z { > 0, i$ E}. 

The interior of the domain is given by d% and the origin is represented by <9{i,2,...,d}- Subsets 
other than the empty set represent the "boundaries" of the state-space and correspond to 
system configurations in which at least one station is empty. The collection of all possible 
values that the increments of the process Q can take depends on the current region at which 
Q is positioned. However, in any case, such collection is a subset of 

V= {ej,-ej + ej,-ej :i,j = 1,2, ...,d}, 

where e* is the vector whose i-ih component is one and the rest are zero. An element of 
the form e* represents an arrival at station i, an element of the form — e* + Cj represents a 
departure from station % that flows to station j and an element of the form — ej represents a 
departure from station % out of the system. The set of all possible departures from station % 

is a subset of 

Vj~ = {w : w = —Ci orw = — e« + tj for some j = 1, ..., d}. 

Because of the nonnegativity constraints on the boundaries of the system we have to be 
careful when specifying the transition dynamics. First we define a sequence of i.i.d. random 
variables {Y (k) : k > 1} so that for each w G V 

F(Y(k) =w) = 



A, 


if w — ei, 


fM*ij 


if w = —ei + ej, 


HiPiO 


if w = — e^ 



The dynamics of the queue-length process admit the random walk type representation given 
by 

Q(k + 1) = Q(k) + C (Q(k), Y(k + 1)) , (8) 

where ( (•) is the constrained mapping and it is defined for x G 8e via 

C (x w) 4 I ° if w G U ^ V * 7 

^ w otherwise 

The large deviations theory associated to Jackson networks is somewhat similar (at least 
in form) to that of random walks. One has to recognize, of course, that the non-smoothness of 
the constrained mapping as a function of the state of the system creates substantial technical 
complications, but we will leave aside this issue in our discussion because our objective is 
simply to describe the form of the necessary large deviations results for our purposes. An 
extremely important role behind the development of large deviations theory for light-tailed 
random walks is played by the logmoment generating function of the increment distribution. 
So, given the similarities suggested by the dynamics of flH]) and those of a simple random 
walk it is not surprising that the logmoment generating function of the increments, namely, 

V(M) = logE [exp {9 T ax,Y(k)))} (9) 

also plays a crucial role in the large deviations behavior of p„ (x) as n /* oo. 

In order to understand the large deviations behavior of p^ it is useful to scale space 
by 1/n, thereby introducing a scaled queue length process {Q n (k) : k > 0} which evolves 
according to 

Q n (k + 1) = Q n {k) + -C (Qn(k), Y(k + 1)) . 

n 

Suppose that Q n (0) = y = x/n and note that 

T {0} = inf {k >l:Q n (k) = 0}, Tl 4 inf {k >l:V(Q n (k)) > 1}. 



Note that using the scaled queue length process one can write 

1 



Pk (y) = E 



P v n (y+-((y,Y(l))) 

n 



(10) 



Large deviations theory dictates that 

Pn (v) = ex P (-nW v (y) + o (n)) (11) 

as n / oo for some non-negative function Wy (•)■ in order to characterize Wy (•) we can 
combine the previous expression together with ( TTOl) and a formal Taylor expansion to obtain 



Pn(y) 



E 



P v n (y + -C(y,Y(i))) 

n 



Eexp{~nW v [y + ~C (y, Y (1))] + nW v (y)} 

Eexp{-dW v (y)({y,Y(l)) + o(l)} 
exp{iP{y,-dW v {y)) + o{l)). 



Sending n/oowe formally arrive at the equation 

ij (y, -dW v (y)) = (12) 

together with the boundary condition Wy (y) = if V (y) > 1. The previous equation 
is the so-called Isaacs equation which characterizes the large deviations behavior of p v (■) 
and it was introduced together with a game theoretic interpretation by Dupuis and Wang 
in [TO]. The solution to (112"]) is understood in a weak sense because the function Wy (■) is 
typically not differentiable everywhere. Nevertheles, it coincides with a certain calculus of 
variations representation which can be obtained out of the local large deviations rate function 
for Jackson networks (see |17j). 

An asymptotic lower bound for Wy (y) can be obtained by finding an appropriate subso- 
lution to the Isaacs equation, in which the equality signs in ( 1121) are appropriately replaced by 
inequalities thereby obtaining a so-called subsolution to the Isaacs equation. In particular, 
Wy (•) is said to be a subsolution to the Isaacs equation if 

V(y, -dW v (y)) < (13) 

subject to Wy (y) < if V (y) > 1. The subsolution property guarantees Wy (y) < Wy (y), 
which translates to an asymptotic logarithmic upper bound p\ (y). The subsolution is said 
to be maximal at zero if Wy (0) = Wy (0). Not surprisingly, subsolutions are easier to con- 
struct than solutions and, as we shall discuss in the next section, beyond their use in the 
development of asymptotic upper bounds they can be applied to the design of efficient sim- 
ulation procedures. The use of subsolutions to the Isaacs equation for the design of efficient 
simulation algorithms was introduced in [10] . A derivation of the subsolution equation ([TBI 
following the same spirit leading to (fT2|) using Lyapunov inequalities is given in [6]. 

As we mentioned in Section 2, the efficiency analysis of a rare-event simulation estimator 
depends on the growth rate of its coefficient of variation. We are interested in an asymptotic 
analysis that goes beyond the error term exp(o(n)) given by the large deviations approxi- 
mation (ITT]) . So, we must enhance the large deviations approximations in order to provide a 
more precise estimate for p\. Developing such an estimate is the aim of the following propo- 
sition which follows as a direct consequence of Proposition 2 and the analysis in Section 5 in 
[I] ; see also Section 4 in this paper for a sketch of the proof. 

Proposition 1. There exists K > (independent of x and n) such that 

p v n (X) < KP{V (Q (oo)) = n}/P{Q (oo) = x}, 

where Qoo is the steady state queue length. Moreover, if \\x\\ < c for some c £ (0, oo) therft 

p v n (x) = n[P{V (Q (oo)) = n}/P{Q (oo) = x}] 

as n / oo. 



2 Given two non-negative functions /(•) and g(-), we say / (n) = £l(g(n)) if there exists c, no G (0, oo) 
such that / (n) > eg (n) all n > n . 
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Remark It is important to keep in mind that we shall mostly work with the process Q (•) 
directly, as opposed to the scaled version Q n (-) which is used in the analysis of [8]. 

The previous proposition provides the necessary means to estimate p% up to a constant; 
we just need to recall that the distribution of Q (<x>) is computable in closed form (see [21] 
p. 95). In particular, we have that 



m jt 



7r(mi,...,m d ) = JJP(Qj(oo) 

3=1 

d 

= n t 1 ~ Pi) p? j > j = - 1 ' •••' d > and m J - °- 



3=1 

We shall use 7r (•) to denote the stationary measure of Q. In simple words, the previous 
equation says that the steady state queue length process has independent components which 
are geometrically distributed. In particular, P (Qj (oo) = m) = pf{l — Pj) for m > 0. The 
next proposition follows directly from standard properties of the geometric distribution (see 
Proposition 3 in [1]). 

Proposition 2. P[V (Q (oo)) = n] — (e~™ 7v 'n /3v '~ 1 ) ; where 7y = — \ogpX , in which p\ = 
max{pj : f j = 1}; and fly = J^ I {pi = pY , v « = 1} is the number of bottleneck stations in the 
target subset corresponding to v. 

4 The Splitting Algorithm 

The previous section discussed some large deviations properties required to guide the con- 
struction of an efficient splitting scheme using the theory developed in the work of Dean and 
Dupuis [8J. In order to explain the construction suggested by Dean and Dupuis let us first 
discuss the general idea behind the splitting algorithm that we shall analyze; a variation of 
which was first applied to Jackson networks by by Villen-Altamirano and Villen-Altamirano 

H5I- 

The strategy is to divide the state-space into a collection of regions {C™ : < j < l n (x)} 
which are nested and that help define "milestone" events that interpolate between the initial 
position of the process and the target set, which corresponds to the region Cq. That is, in 
our setting we put Cq = {y G S : V (y) > n} and the remaining C"'s are placed so that 
Cq C C™ C ... C Cj^ n . How to construct the level sets C- in order to induce efficiency will be 
discussed below. An observation that is intuitive at this point, however, is that one should 
have M n = (n) so that the next milestone event becomes accessible given the current level. 
For the moment, let us assume that the CJ's have been placed. The splitting algorithm 
proceeds as follows. 

Algorithm SA 



3 Given two positive functions / (•) and g (•), recall that / (n) = 9 (g (n)) if / (ft) = O (g (ft)) and / (ft) = 
n(g(n)). 
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1.- Initiate the simulation procedure with a single particle starting from position x G C^ 
for a given k > 1. Let u»i = 1 be the initial weight associated to such particle. 

2.- Evolve the initial particle until either it hits {0} or it hits level C^_ v If the particle 
hits {0}, then the particle is said to die. If the particle reaches level C^_ x then it is replaced 
by r identical particles (for a given integer r > 1). The replacing particles are called the 
immediate descendants or children of the initial particle, which in turn is said to be their 
parent. The children are positioned precisely at the place where the parent particle reached 
level C^_ x . The weight Wj associated to the j-th children (enumerate the children arbitrarily) 
has a value equal to the weight of the parent particle multipled by 1/r. 

3.- The procedure starting from step 1 is replicated for each of the offspring particles in 
place; carrying over the value of each of the weights at each level for the surviving particles 
(the weights of the particles that die can be disregarded). 

4.- Steps 1 to 3 are repeated until all the particles either die or reach level Cq. 

Dean and Dupuis in [8] show how to apply large deviations theory to select the C^'s in 
order to obtain a weakly efficient splitting algorithm. One needs to balance the number of 
the C^'s so that it is not unlikely for a given particle to reach the next level while keeping 
the total number of particles controlled. We now provide a formal motivation for the use of 
large deviations for constructing the C"'s in a balanced way. 

It is convenient, as we did in our formal large deviations discussion in the previous section, 
to consider the scaled process Q n (•). Let us assume that the splitting mechanism indicated 
in Algorithm SA is in place and that our initial position is set at level Q (0) = x, so that 
Qn (0) = y — x/n. The C™'s are typically constructed in terms of the level sets of a so-called 
importance function which we shall denote by U {■). In particular, put D n = {y £ n~ l S : 
V (y) < 1} and set C™ = nL Zn ^, where 

L z ±{yeD n :U(y)<z}, (14) 

and the z n (j)'s are appropriately chosen momentarily. Then, define 

l n (y) = min{j > : ny e Cf} = min{j > : x E Cf}. (15) 

The total weight corresponding to a particle that reaches level Cq given that it started at 
level l n (y) is r~ ln( ~ y \ In order to have at least a weakly efficient algorithm we wish to achieve 
two constraints. The first one imposes the total weight of a particle reaching level Cq to be 
p^ (x) exp (o (n)); this would guarantee that the second moment of the resulting estimator 
achieves asymptotic optimality. The second constraint dictates that the expected number of 
particles that make it to Cq, which is roughly r ln ^p^ (x) exhibits subexponential growth (i.e. 
exp(o(n))); this would guarantee a cost per replication that is subexponential. Note that 
both constraints lead to the requirement of r ln ^p^ (x) = exp (o (n)). So, given a subsolution 
Wy (■) to the corresponding Isaacs equation, which implies that 

Pn ( x ) < ex P {—nWy (y) + o(n)) , 

it suffices to ensure that 

l n (y) log (r)-nW v (y) = o(n). (16) 
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The behavior of l n (y) as a n /* oo only relates to the properties of the function U (•) and it 
is really independent of the large deviations behavior of the system. In particular, picking 
z n (j) = Aj/n, A G (0,1] yields Z n ($/) = [rat/ (?/) /A] and therefore, equation ([16]) suggests 
that one should select U (y) = AWy (y) / log (r) with Wy (0) = Wy (0) in order to obtain 
a weakly efficient estimator for p^. This is precisely the conclusion obtained in the work 
of [8] who present a rigorous analysis that justifies the previous heuristic discussion. Our 
development in the next section will sharpen the efficiency properties of the sampler proposed 
in [8] when applied to Jackson networks. So, we content ourselves with the previous heuristic 
motivation for the splitting method that we will analyze in the next section and which in 
turn is based on the viscosity subsolution given by 

W v (y) = Q T y-logpV, (17) 

where Qi = — logpj for i — 1, ..., d, see e.g., [TT] and [8]. 

We close this section with a precise definition of the estimator that we will analyze. First, 
given a constant AG (0, 1] (the level size) define Wy (•) as indicated in (ITT)) for each y = x/n 
with x G S. Then, select an integer r > 1 and define U (y) = AWy (y) /log (r). Given the 
initial position x put y = x/n and define the sets {C™ : 1 < j < l n (y)} as indicated above 
(see equation (1151) ). Run Algorithm SA and let N n be the number of particles that survive 
up to Cq; their corresponding final weight is l/r ln( - x \ Our estimator for p^ (x) is simply 

R n (x) = N n (x) /r l ^ x) (18) 

Now, for the sake of analytical convenience, when analyzing the second moment of R n (x) 
we will adopt the so-called fully branching representation of the previous estimator (see [8]). 
Such fully branching representation is obtained by splitting death particles at level zero. It 
is useful to think about fully branching conceptually in the following terms. Start with a 
particle at position x in level C" , -, and run step 2 of Algorithm SA, but instead of killing the 
particle if it hits {0} before reaching C^,^^ just allow it to also split into r particles and 
update the weight of the children as indicated in Algorithm SA when the particle reaches 
^? n (x)-v Step 3 continues, now the particles that are sitting in level C] l , x - ) _ 1 will evolve and 
the death particles will once again split and remain in state (so state is being populated 
with death particles). After l n (x) iterations we have r ln( - x > total particles labeled 1, 2, ...,r ln( - x \ 
each with weight l/r ln( - x '. We define Ij as the indicator function of the event that the j-th 

particle is in Cq so that N n (x) = Y^j=i Ij- ^ ne ^ u ^y branching representation of R n (x) is 
simply 

|Jn(«) 

Rn (x) = r- 1 ^ J2 l r ( 19 ) 

5 Analysis of Splitting Estimators 

We are now in a good position to perform a refined efficiency analysis for the estimator R n (x). 
We shall break our analysis into two parts. The first part corresponds to the expected number 
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of particles generated per run and the second part deals with the second moment of R n (x). 
First, we will review a technique studied in [1] based on the corresponding time reversed 
chain associated to the queue length process Q. For both quantities we are able to obtain 
an upper bound. We are then able to reach the conclusion that this multilevel Splitting 
algorithm substantially outperforms the direct polynomial time algorithm for solving the 
associated system of linear equations. 

Our analysis takes advantage of the time reversed process associated to the underlying 
Jackson network which we shall now define. Given the transition matrix {K (x, y) : x,y G S} 
of the process Q, we define the reversed Markov chain Q = {Q (k) : k > 0} via the transition 
matrix K (•): 

K (y, x)=K (x, y) vr (x) /it (y) , 

for x,y G S. It turns out that Q also describes the queue length process of an open stable 
Jackson network with stationary distribution equal to vr(-), (see [21] p. 95). We will use 
P x (•) to denote the probability measure in path space associated to Q given that Q (0) = x. 
The following result is similar to that of Proposition 2 in [1] . However, our representation 
in (|20|) is slightly more useful for our purposes. 



Proposition 3. 



v K (Q (0) G Cft, f {x} < f {0} , f {x} < f c ^ 

7r(x)P x {T {x} > T c % A T {0 }) 

*\ 

(21) 



p f Q(o)6Cj,r w <r { o } <Tcj 



71 



x)P (T {x} <Tc S AT {0} ] 



where f cs = ini{k > 1 : Q (k) G C£}, f {x} = mf{k > : Q (k) = x}, T cs = ini{k > 1 : 
Q (k) G Cq} and T{ x y = inf{k > : Q (k) = x}. Moreover, there exists 5 > (independent 
of x and n) such that 

P X {T {X} > T cs A T {0} ) > P X {T M > T {0} ) > 5, (22) 

where \\x\\ is the L\ norm of x. 

Proof. We assume that i^O. The case x = is included in the analysis of ( 12"T|) . First, we 
observe that 

Pn (?) = P x (Tcj < T {0} , T {x} < T cs A T {0 }) + P x (T cs < T {0} , T {x} > T cs A T {0} ) 
= p v n (x) P x (T {x} < T cs A T {0} ) + P x (T cs < T {0} ,T {X} > T C n A T {0} ) . 



Therefore, 



ry m - - x ( Tc o <; T W' T w - Tc s A T i°}) 

Pn \ x ) ~ 



Px (T{ x } > T C n A T{ }) 
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Following the same technique as in Proposition 2 in [I] we have that 
vr (x) P x {T C n < T {0} , T w > T cs A T {0} ) 

CO 

= ^tt(x)P x (T C n < T {0} , T w > T C n A T {0} , Tq™ = k) 



k=0 
oo 



^tt(x) ^ K (Vo,yi) x - x K (Vk-i,Vk) I (k < T {0} ,T {x} > k AT {0} ,T CS = k) 

fc=i yo=x,yi,~,Vk 



OO 



= J2 J2 K ' (yi> y°) x - xK ' (y*> ^-0 * fo*) J (^ < r {0} , t w > ^ a t {0} , t cs = k) . 

fe=l y =x,y 1: ..,y k 

Letting ?/,• = ?/fc_j for i = 1, ...k we see that the summation in each of the terms above ranges 
over paths y' , ■■■,y' k satisfying that y' G Cq, T{ x y = k (so in particular y' k = x) and also that 
7{o} > k, Tc™ > k. So, we can interpret the previous sum as 



K (q (0) e C n , f w < f {0} , f w < fcs) ■ 



This yields part f[2"Uj) . Part ( 12~T|) corresponds to Proposition 2 of [1] . Finally, ( 12 2 j) also follows 
as in Proposition 7 of |1]. D 

Proposition 1 and 2 from Section 2 follow as a consequence of this result, the rest of the 
details are given in Section 5 of [4J. 

Given the subsolution we proposed in Section 4, the importance function can be written 

as 

U (x/n) = Wy (x/n) = I — q t x — logpjf 



log r \n * J log r 

= C f A - -a T xA J , (23) 

where C = — log p^/ log r, and a = g/logpY- The level index function also simplifies to 

r (x) = [ n[/ ^ /n) l = fnC (l - -a T x)^ = \C(n- a T x)]. (24) 



We shall first look at the expected number of surviving particles of the splitting algorithm 
which characterizes the stability of the algorithm. One shall keep in mind that when the 
complexity of the splitting algorithm is concerned, what actually matters is the total function 
evaluation involved in each run. An upper bound is obtained for this quantity, as measured 
by the sum of all particles generated at interim levels weighted by the maximum remaining 
function evaluations associated with each of them. We first have the following result. 

Proposition 4. The expected terminal number of particles for the splitting algorithm specified 
by (A, U) above satisfies 

E[N n (x)] = e(n^ v ~ l ) (25) 

where /3y, introduced in Proposition 2, denotes the number of bottleneck stations correspond- 
ing to the vector v. 



15 



Proof. Note that 



E[N n (x)] = r in WpV( x ). 



In Proposition [2] we know that p. 



v 



X 



-C 



(e-T" V^ -1 ) . Since e" 



IV 



, we can write p% (x) = (r nC n l3v *) . Hence, plug in / n (x) = \C C 
E [N n (x)} = (r-nC n p v -i r \nCJ}j = q ( n /V-i) _ 



e 

n - 



T 



-Clogr 



) ] , we have 



D 



As pointed out earlier, the number of terminal surviving particles, although serves as a 
reasonable proxy to measure the stability of the algorithm, is not suitable for quantifying the 
complexity We also need to take into account the number of function evaluations required 
to generate R n (x). The next result addresses precisely this issue. 

Proposition 5. The expected computational effort per run required to generate a single 
replication of R n (x) is O (n^ v+l ) . 



Proof. To see this, let iVT' 



/;/ 



0, ...., l n (x), be the number of particles that survive to level 



rin 



(x)— to' 



Also let r] m j be the remaining computational effort of the j-th particle at the start 



of the m-th level until it either reaches the next level or it dies out. Put r] m j (xj) to be 
the expectation of r\ m j given that the position of the j-th particle at the start of level m is 



x 



j- 



Note that the norm of the position of Xj is less than c • m for a given constant c that 
depends on the traffic intensities of the system but not on the position of the particle per-se. 
Therefore, it is easy to see that 

sup fj m ,j (xj) <c-m, 

for some c G (0, oo). Intuitively, each particle at level m either advances to the next level, or 
it dies out by hitting the zero level before moving to the next one, since it takes © (1) work 
to cross one single layer, r) m> j is dominated by the work required to die out, and hence its 
mean is bounded from above bycxm for some constant c. The expected total work per run 
is then given by 



E 



l n (x)-l AT™ 
m=0 j=l 


l n (x)-l 

= E E 

m=0 


r N n i 

/ j Vm,j \ x j ) 






ln{x)-l 

< J2 Era-cm 

m=0 




ln(x)-l 

m=0 




l n (x)-l 

*°- E 


®* = o 


(n^ +1 ) , 



m=0 



for some positive constant c where the first inequality is due to independence. 



□ 
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To facilitate the analysis of the second moment of R n (x) we add the following notations. 
We follow the analysis in [8j to make our exposition here self-contained. For a given generation 
m, denote by Q m j the position of the jth particle; recall that the accumulated weight up 
to the mth stage of such a particle is r m . Let Xmj be the disjoint grouping of particles in 
the next generation (i.e., m + 1) according to their "parents" in generation m. For k e Xm,j, 
denote by dk the offsprings of this particle at the final stage l n (x). We then have the following 
expansion of the second moment of R n (x): 



E, 




ln(x)-l 

E E - 



m=0 



E E Ew 

j=\ k,l£Xm,j,k^l \m k £d k 



-ln(x) 



(26) 



J2 Im ' 1 



-ln(x) 



+ E X 



E'> 



-2l n (x) 



where we define I mk to be the indicator function of the event that particle m^ is in the set 
Cq. The second term above is essentially the diagonal terms of the second moment ([2"B"j) . 
and for the off-diagonal terms, for each generation, we categorize particles according to their 
common ancestors, a technique used by [5|. For the first term, we have 



l n (x)-l 
m=0 

ln(x) 

m=0 



E E Ew 

- r m 



-ln{x) 



E* 

\miGdi 



-ln(x) 



E (7EW 



-(ln(x)-m-l) 



fc)'GXm,.7 i^T 1 '' 



m k £d k 



;E* 

rrnedi 



-{l n {x)-m-l) 
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Conditioning on the whole genealogy up to step m, we obtain 
J2HV(Qm, 3 )>0){r- m ) 2 

E (\ E W-«wi— ") (i £ /, 

J2nV(Q m ,)>0){r- m ) 2 



—(l n (x)—m—l) 



i=i 



E ( E f; E w-<m.>— .>] (1 53 /„,, 

\k,l€Xm,j,k^l \ m k £d k / \ mi&di 



-(ln.(x)—m—i) 



1L 



^/(l/(Q mJ )>0) 



-2m 



E (;*** ( E ^r-<^- m -A Ie^ ( e '- 

Xm,j,k^l \ \m k ed k / \mi£di 



—(l n (x)—m—l) 



k,l&Xm,j,k=/=l 



Note that E Qm> . [E mfc&Ifc / mfc r-('»(*)-m-D] = p£(Q«y). and W = £^r~ 2 = r/(l-r). 
Adding over m we obtain 

l n (x)-l r r m 

m=0 Li=i 

/ w (a)-l 

= Wj] r- m E x [I {V (Q m>1 ) > 0)^ (Q mj ) 2 ] . 

m=0 

Combining this with the diagonal term in (126|) . which can be readily expressed as r~ in "p^ ( 



,r 



we arrive at the following expansion for the second moment of R n (x) 

ln(x)-l 



E x [R n (x) 2 } = W J2 r- m E x [l(V(Q mA )>0)pV(Q mA ) 2 ]+r ->■',,;,[.,■). 



,-l n (x) V 



(27) 



m=0 



The next result takes advantage of expression ( !27l) to obtain an upper bound for E x [R n (x) ] . 
Proposition 6. The second moment of R n (x) satisfies 

E [R n (x)} 2 = v v n {xf O (n?) . (28) 

where j3 = Ei=i 1 {pi = P*) is the number of bottleneck stations in the whole network. 



In order to prove the previous result, we will show that the second moment of R n (x) 
is dominated by the first item on the righthand side of the equality in (1271) . In turn, the 
asymptotic behaviour of such term hinges on the conditional distribution of the exact position 
of the particle in generation m, Q m ,i in CJl ( x )- m - 

Proof. We begin the proof with an important property implied by the splitting algorithm: 

V (Q m ,i) > <^> Q m>1 G Cl b n(x) _ m = nL( ln{x) _ m)A/n 

&■ Q m ,i E {z E nD n : U (z/n) < (l n (x) — m) A/n} 

& Qm,i e {z G nD n : C I A a T zA J < — (\C {n - a T x)] - rn)} 

=>- Qm i G {z G nD n : C ( 1 a T z ) < - (C (n - a T x) - m + l) } 

\ n J n v v 



^ Qm,i G {z G nD n : a T z > a T x H — — } 



C 



& Qm,i ^ {z G nD n : 2 < £ x — (m — 1) logr} 



(29) 



where we used the representations of U (•) and Z n (x) in (|23|) and (1241) and the definition of 
L 2 in ( Tl4l) . In other words, if a particle survives m generations then its current position 
is beyond the mth level, which implies that the weighted sum of system population, with 
weight given by the vector g, is bounded from above by that of the initial position adjusted 
by a linear function in m. If we define the stopping time Tm = inf{k > 1 : a T Q (k) > 
a T x + in ^-} = inf{k > 1 : g T Q (k) < g T x — (m — 1) logr}, the above property also implies 
that V (Q m ,i) > => Tm. < T . Now, the expectation term in the sum of ( l2"Tj) can be 
expressed as 



^ X [l{V{Q m , l )>Q)p V n {Q m ,lf] 

<E X [I (g T Q m>1 < g T x - (m - 1) logr) p v n (Q mA ) 2 ] 

=E X [p v n {Q m>1 f \g T Q m>1 < g T x - (m - 1) logr] P x (f f < T 



=E 3 



P^ T v n < T \Q fm \T~ < T , £ 7 Q m>1 < ^ x - (m - 1) logr 



P, ( T™ < T 



(30) 



where we used the property derived in fl29|) . For the first item in (130]) . we have 



E, 



F 2 X T v n < T \Q fm |r™ < T , ^ Q mil < ^ x - (m - 1) logr 



<KE 



vr 2 (Cg 
Ltt 2 (Q. 



|^ T Qm,i < Q T x - (m - 1) logr 



m.l. 



,ni2. 



<ci K v_1 BT E » e_2e gm,1 l0 T Qm,i < T * - (m- l)logr 



(31) 



where ci, i^ are some constants independent of n. Here we used Propositions 1 and 2 for the 
the last two inequalities respectively. As for the expectation term in (13"T1) . since the process 



19 



Q (•) has for each dimension an increment at most of unit size, we can write 



E w 



-28 1 



,1 \Q T Qm,i < Q T x - (to - 1) logr 



E„ 



-2 e 3 



1,1 1 q x — (to — 1) logr — 5 < £> Q m ,i < £ £ — (to — 1) logr 



< C2 exp (— 2^> x + 2 (to — 1) log r) 



c 3 exp 



2 m Jn ^^"\ -„_( n v\- 2S rr 



C 



logplf =c 3 (pD 



(32) 



where 02,03 and 5 are some positive constants. It remains to give a bound on the term 



v T™ < T ) . As a result of Proposition 2, 



P* ( % < T 



1 

7T (x) 
1 



< — ^rP [e T Q (oo) < g T x - (to - 1) logr] 

-P 



7T X] 



d T Q(oo) > d T x + ^" " ' 



C 



where C = - log p*/ logr, p* = max^pj e (0, 1) and d = g/p* = (logpi/logp*, ..., logp rf /logp* 
Note that dj G (0, 1). To finish the proof we need the following Lemma. 



Lemma 1. 



P 



d T Q (oo) > a T x + 



(to — 1) 



e 
e 



p Z(/?,1-P*)> 



TO — 1 



TO — 1 



3-1 



(P* 



m — 1 

c 



where Z (n,p) denotes a NBin (n,p) (negative binomial) random variable. 
Proof. Note that 



d Q (oo) = Q (oo 



,t e 



logp* 



^ <2i (oo) / (pi = p*) + ^ Qi (oo) / (pi 7^ p*; 



i=i 



log Pi 
logp* 



= Z(/3,l-p,) + W- 
One direction is elementary, since a T Q (oo) > Z (ft, 1 — p*), we clearly have 



P 



d T Q (oo) > a T x 



(to — 1) 



> P 



Z (f3, 1 — p*) > d T x 



(to — 1) 



(33) 
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For the other direction, note that there exists a constant p < p* such that 



w = J2QiMHPi^ 



i=i 



log Pi 
logp* 



< J2 Qi (°°) l iPi 7^P*) = Z(d-(3,l-p). 



i=i 



As a result, 



A.T 



a 1 Q (oo) < st Z((3,l-p*) + Z(d-[3,l- p) 



where " < st " denotes that the left hand side is stochastically dominated by the right hand 
side. But since 1 — p* < 1 — p, a similar argument as given by Proposition 3 in [4] allows us 
to obtain 



P 



a T Q oo > a T x + - — s — - 

C 



<c P 



Z (/3, 1 — p*)> a. x 



T (m—l) 



C 



(34) 



for some finite constant Cq that is independent of m. Combining (I3"3"j) and (13~4"|) . we have 



P 



~t^ / x -t (m — 1) 



6 



P \Z(P,l-p*) >a 1 x + 



r V m — 1) 



c 



(35) 



Using again Proposition 3 of [4], we reach the conclusion that 

(m — 1 



P 



r. 



d Q (oo) > a x + 



C 



6 



m — 1 \ 



C 



(p*y 



n 



Going back to (130]) . the above Lemma allows us to write 



W> x [Trn<T ) <C 4 



m — 1 \ i 



C 



(p*y 



(36) 



If we combine ( |3TT) . (!32l) and ( J36l) . we obtain the following upper bound for the expectation 
term in the sum of expression ( T27j ) : 



E, [/ (V (g m;1 ) > 0) p£ (Qm,i) ] < cj£ (*) MO 



1 1 m — i 



c 



p-i 



(p* 



m—l 

c 



CP n ( X ) {P* ) 



c 



0-1 



V / \2 rn-1 



CPn ( X ) r 



m—l 

C 



0-1 



(37) 



where for the first equality we use the fact that pY < p* and C < C, and for the second 
equality we use pX = r~ c . Putting the bound in ( 137|) back to the sum in the first item of 
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271). we have 



*n(x)-l 

^ r- m E x [/(V( mil )>0)^( m)1 ) 2 ] 

m=0 

'n(x)-l / 1 \ /3-1 



< cr" 1 ^ pi (xf I ^-f— ) 

m=0 V C- / 

= p£ (^) 2 O (n?) (3* 



Finally, note that the second item of (12"T|) is dominated by ( 13"8|) , and it follows immediately 
that 

E[R n (x)] 2 =p^(x) 2 0(n p ). 

D 

Equipped with these results, we are ready to summarize our discussions in the statement 
of the following Theorem, which is the main result of this paper. 

Theorem 1. To estimate the overflow probability p^ (x) using R n (x), the number of function 
evaluations needed for a given level of relative error is O (n 1 



Pv+P+A _ 



Proof. Recall from section [2] that the number of function evaluations sufficient to achieve a 
pre-determined level of relative accuracy for the Splitting estimator is proportional to the 
work-normalized squared coefficient of variation. This is therefore immediate by combining 
the upper bound analysis of the computational effort per run in Proposition [5] along with 
the upper bound of the second moment of R n (x) available in Proposition [6j □ 

A direct comparison to the O (n 3d ~ 2 ) complexity of solving a system of linear equations 
(see Section [2] ) yields the immediate conclusion that the Splitting algorithm is "efficient" in 
the sense that it improves over the "benchmark" polynomial algorithm. Even in the worst 
case when we look at the total population of the network and the network is totally symmet- 
ric, i.e., all stations are bottlenecks (j3y = P = d > 3), the number of flops needed is off by a 
substantial factor n d ~ 3 . In the case where /3y = fi — 1, the algorithm only requires a number 
of function evaluations that at most grows cubically in the level of overflow n. Furthermore, 
if the number of bottlenecks is less than half of the total number of stations, i.e. j3 < d/2, the 
Splitting algorithm enjoys a running time of order smaller than O (n d ) , which is not worst 
than storing the vector that encodes the solution to the associated linear system. If, on the 
other hand, more than half of the stations are bottlenecks, faster importance sampling based 
algorithms do exist at least for the case of tandem networks; see the analysis in [6], which 
implies that O (n 2 ^~^ +1 ) function evaluations suffice to obtain an estimator with a given 
relative precision. Overall, the analysis thus provides some sort of guidance on the choice 
of simulation algorithms. It is meaningful to point out that the previous comparison is not 
based on the sharpest analysis. In fact we only resort to a rather crude upper bound in the 
analysis of the second moment of R n (x) in (I3~T1) . A sharper result is possible by bounding 
the expectation term in (1201) with more care. But as pointed out in the Introduction, even 
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though there is still room for a more refined analysis, we believe our work provides substan- 
tial insights leading to a better understanding of the relations between these two classes of 
algorithms. 
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